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' Leaves and flowers frequently have a characteristic rippling pattern at their edges. Recent experiments found 

, similar patterns in torn plastic. These patterns can be reproduced by imposing metrics upon thin sheets. The 

■ goal of this paper is to discuss a collection of analytical and numerical results for the shape of a sheet with a 

' non-flat metric. First, a simple condition is found to determine when a stretched sheet folded into a cylinder 

loses axial symmetry, and buckles like a flower. General expressions are next found for the energy of stretched 
sheet, both in forms suitable for numerical investigation, and for analytical studies in the continuum. The bulk 
) of the paper focuses upon long thin strips of material with a linear gradient in metric. In some special cases, the 

^— ( I energy-minimizing shapes of such strips can be determined analytically. Euler-Lagrange equations are found 

. which determine the shapes in general. The paper closes with numerical investigations of these equations. 

<N : 

PACS numbers: 45.70.Qj,02.40.-k 
O ' I. INTRODUCTION 

c/3 ; 

d , A characteristic rippled pattern often appears at the edges of leaves and flowers. One can even produce it by ripping in half a 
^ ' thin sheet of plastic, such as a plastic bag[l|. A numerical example of a stretched sheet with two generations of waves appears 
I in Figure |l]. I have found a case where the mathematics of this problem can be solved exactly. 
■ O ■ This topic touches on an interesting basic question, which concerns the complexity of the instructions needed to generate 
^ ' complex natural patterns. The dominant belief among biologists is that the curling shapes of plants are produced by detailed 
genetic instructions telling various sections to curl up and down|Q|. The physics community interested in studying patterns 
would prefer to assume that complex forms are produced when possible by simple rules, a view with some support also among 
biologists ||^ ^. It is natural to guess that by imposing non-flat metrics upon thin sheets, elasticity alone will compel them to 
curl spontaneously into fractal forms. However, to proceed from a guess to a detailed demonstration requires some effort. 

The basic question addressed throughout this paper is the following: Suppose one has a rectangular sheet of material, described 
by coordinates x and y. Take this rectangular strip and impose a new metric on it so that distances dr between nearby points 
originally separated by {dx, dy) are given by 

OO ■ 

O' dr^ ^g^,dx^ ^-Ig^ydxdy + gyydy^. (I) 

(N ' 

, Embed the sheet in three-dimensional space, allowing it to curl as needed to obey this new metric. What shape does it take? 
' The main results in the paper are the following: 

a \ 

^ \ • Use of results from differential geometry to find when axisymmetric sheets wrapped into cylinders are stable against 

I ■ buckling (Section II). 

; 

' • Energy functional for a model of a thin sheet with non-flat metric as a collection of interacting mass points (Section III). 

o : 

O ■ • Derivation of continuum nonlinear elastic theory from the original energy functional (Section IV) 

. ^ ' • Definition of solvable problem in thin strip with Unear gradient in metric in terms of a set of ordinary differential equations 

^ ! (Section V). 
' 

' • Exact solution of these ordinary differential equations in terms of elementary function when certain bending and torsion 
angles in them are specified in specific ways (Sections VI-VIII). 

• Calculation of continuum energy functional in terms of the solutions of the ordinary differential equations (Section IX). 

• Derivation of Euler-Lagrange equations from continuum energy functional, demonstration that analytical solutions of 
previous sections coiTespond to energy minimizing solutions when special boundary conditions are imposed (Section X). 

• Numerical solutions of Euler-Lagrange equations to find shapes of strips for more general boundary conditions (Section 
XI). 
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FIG. 1 : Illustration of shapes produced by stretching a membrane at one end. This figure was produced by numerical minimization of Eq. [lo| 
with gxx (y) = .7exp(— j/12) + .3exp(— j/36). The numerical sample was 600 units long, 60 units wide, and two layers high. 



Nechaev and Voituriez[H] and Henderson and Taimina |^] have carried out the only other studies along these lines of which I 
am aware. The first group of authors provides an exact solution using conformal mapping techniques to analyze an exponentially 
growing metric, while the second set of authors provides additional information on the mathematics of hyperbolic planes. Both 
studies are nicely complementary to the work presented here. 

II. BUCKLING OF FLOWERS 

It would be very appealing if one could employ differential geometry alone to resolve questions about the buckling of sheets 
with non-flat metrics. I have not had success along these lines for the problem depicted in Fig. |l]. However, if I digress briefly 
to discuss flowers rather than leaves, I can make some progress, and want to mention this result before passing on to the main 
topic of the paper. 

Suppose that one adds an additional constraint to the problem, which is that when the sheet is embedded in three dimensions, 
it is required to be periodic along x, with period InRoo- One has a distorted cylinder, as shown in Figs. ^. One can use results 
from differential geometry to obtain a simple criterion for when such a structure loses its cylindrical symmetry and buckles. The 
Gauss-Bonnet theorem[p| applied to this case says that 

dsK = - I dAK, (2) 



where k, is the geodesic curvature, and K is the Gaussian curvature. The first integral is a line integral taken around the edge of 
the cylinder, while the second is a surface integral. 

Use X and y to refer to material coordinates. For a metric where only gxx{y) differs from the value expected in flat space, the 
Gauss-Codazzi relations!^ give 

^=-^'-W- (3) 

Suppose that after embedding the sheet retains cylindrical symmetry. Then along the outer rim of the flower, k is constant. 
One can place a bound on n. It cannot be greater in absolute value than \/R, where R is the radius of curvature at the rim, and 
it obtains this value only when the edge of the cylinder has splayed out so that the whole outer boundary shares a single tangent 
plane containing the boundary. Thus the left hand side of Eq. || obeys 

-27r< J dsK = 2nRK<2TT (4) 

Turning now to the right hand side of Eq. ^ one can perform the integral in cylindrical coordinates. At any point indexed by 
y, the radius of the cylinder is Rooy/gxx, where Roc is the radius of curvature at the undistorted end of the flower. Along the y 



FIG. 2: Result of numerical minimization of Eq. [lO| for system twisted around into a cylinder The x coordinate of the original material travels 
in an angular direction around the axis, while the y coordinate of the original material describes motion along the axis, (a) The metric is 
y/gxx ~ 1/(1 + 27ry/300), and the original radius of the cylinder is ^ = 47.8 = 300/(27r). The product /?oo«: in Eq.|^ equals 1. The edge is 
without ripples and has splayed out to lie in a plane perpendicular to the axis, (b) Same as (a) except that the metric is ^gxx = 1/(1 +y/2). 
The bound in Eq. ^ is exceeded by a factor of 25. There are around 22 ripples at the edge. 

direction, since gyy = 1, distances along the sheet are unchanged by introduction of the metric, and one can write the integral 
over the Gaussian curvature as 




(5) 



(6) 



= ~2tiR, 



OC 



dy 



Returning to Eq. ^ one has the following condition: 



1^. 



■OO 



dy 



I<1- 



(8) 
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For metrics with slopes less than this value, it is possible for a cylindrical sheet to maintain cylindrical symmetry, but when the 
gradient of the metric becomes too steep, a cylindrical sheet must begin to buckle. 

Fig. ^a) shows such a sheet at the point where the metric has been chosen to make the inequality in Eq. ^ into an equality, 
and the edges of the cylinder are splayed out as far as they can go without buckling. Fig. ^(b) provides an example of a buckled 
cylinder, and the result suggests that the number of ripples at the edge of a cylinder might be given roughly by the quantity on 
the left hand side of Eq. ^ I have not pursued this idea further. 

III. ENERGY OF A THIN SHEET 

I will now proceed to consider leaves, which is to say rectangular strips of stretched material. I will present some particular 
cases where the shapes of rectangular strips with non-flat metrics can be determined, even including very special cases where 
the shape is described in terms of elementary functions. To begin, I will establish in general the energy of a thin sheet with a 
metric that is not flat. 

The energy of a thin sheet is conventionally given by the Foppl-von Karman equations, and is the sum of two terms, one 
involving bending of the sheet and the other involving stretching [||^]. This starting point is inconvenient for three reasons. First, 
many formulations employ coordinate systems that are not general enough to encompass the folds and overhangs that occur in 
this problem. Second, I will deal with sheets that have permanently been stretched, and the equations must be generalized to 
encompass the deformation. Finally, I need to move easily back and forth between numerical and analytical approaches, and the 
numerical discretization of the conventional equations is not simple. 

I avoid all these problems by starting with a physical model of a sheet based upon a discrete collection of interacting points. 
I derive the continuum theory from the discrete model rather than by discretizing continuum equations. Similar numerical 
techniques have been employed frequently in studies of crumpled paper and tethered membranes |llC|, |ll|] . 

Let M, be a collection of mass points that interact with neighbors, and at rest form a thin flat sheet. Let A,j be equilibrium 
vector displacements between neighbors / and j. When the neighbors are not in equilibrium, the distance between them is 
Uij — \uj — M, |. Take the energy corresponding to locations of the mass points to be 

^ = ^E[4-A?;f, (9) 
('7> 

where /C has dimensions of energy per volume, a has dimensions of length, and the sum is over pairs of neighbors. 

With Eq. ^ as a starting point, it is extremely easy to see how to modify equations of elasticity to incorporate a new metric. 
Write 

^=^Y>^j-T.^S.,Af^, (10) 

where g^p is the metric tensor describing deformations of the sheet. All that has happened, in short, is that the equilibrium 
distance between material points has changed. My hypothesis is that by specifying different metric tensors g, often ones with 
very simple functional forms, one can describe all the buckling cascades we have observed in the laboratory. In particular, since 
the plastic sheets are torn uniformly in the x direction, I will assume that the metrics g are constant in the x direction, and vary 
only along y. Often I will take g to be diagonal, with gxx{y) a smoothly varying function and gyy — I. 

IV. CONTINUUM LIMIT 

The long-wavelength deformations of the points m, correspond to elastic deformations of a continuous sheet. To form this 
correspondence, make the replacement 

uj^u{n) + {Aij-V)u{n), (11) 

where r, is the location of m, in an equilibrium of Eq. |[ Define the strain tensor 

ea/3 = 5 



du"^ du'' 
^ dra drp 



(12) 



This definition reduces to the strain tensor of linear elasticity for small deformations and flat metrics, and generalizes it appro- 
priately when deformations are large and the metric tensor g differs from the identity. In terms of the strain tensor, £ can be 
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rewritten as 



/ J nbr. of ( 



E 

a/3 



(13) 



One has to specify a particular lattice in order to proceed further. If one takes the points m, to sit on two two-dimensional 
triangular lattices, stacked over one another as in the first stage of forming an hep lattice of lattice constant a, the the energy 
takes the particular form 



19[e„ + e,v]' + 38[el, + e^, + 2e2,,] 
-32e^j+16[evv + ejez, 
_ + 32 [e% + tl} + 8 \^[e„ (evv - e,,) + 2t„t,,] 



(14) 



I will employ this continuum functional later to develop analytical criteria for membrane shapes. Not all terms in it are equally 
important. However, before arriving at specific deformations to insert into this functional, it is difficult to tell which terms are 
large and which are small. I therefore begin with a geometrical description of wrinkled sheets, and then return to the question of 
which shapes minimize energy. 



V. SOLVABLE PROBLEM AT EDGE OF STRIP 



The full problem at hand is find the minimum energy state of Eq. 10 or Eq. |14jfor a very long strip of finite width and very 
small thickness f, subject to a metric gxx{y) where g„ has some value go aXy ~Q (left hand side of lower panel in Figure [l]), and 
decreases monotonically toward 1 as y approaches the other side of the strip (right hand side of lower panel in Figure nl). It seems 
very unlikely that this problem has in general an exact analytical solution. It is not even clear at the outset whether minimum 
energy states inserted into the functional Eq. |l^ should produce the energies proportional to t that would be characteristic of 
stretching, or energies proportional to that would be characteristic of bending. 

Therefore, it is useful to find cases where the problem can be solved exactly. The motivation for the solvable problem comes 
by looking at the lower panel in Figure |l} Imagine making a new infinitely long strip by slicing off the material a short distance 
w from the left hand side. If w is small enough, then the metric g„ within it should have the form of constant go plus a term 
linear in y. If the original length of this strip along the x direction was L, now its arc length is y/goL. However, in order to be 
able to join onto the rest of the strip on the right hand side, its total extension in the x direction must remain L. 

Therefore, I study a thin strip whose metric is 



V8^ = V8^ii-y/R); fe. = i (15) 

where ^ is a constant. I look for solutions subject u{x,y,z) subject to the constraint that there be some period A for which 

u{x + \,y,z) = xX + u{x,y,z). (16) 

To obtain the benefits of symmetry, also take 

ye[-w/2,w/2] (17) 

so that the center line of the strip is at y = 0, rather than one edge. 

If a long strip of material is given metric Eq. |l5] and no constraint is applied, then the minimum energy configuration is easy 
to find. The material can relax completely by forming a ring of radius R, which curls round and round in a circle. Thus one can 
obtain some intuition about the problem by cutting out the paper figure in Figure || and pulling the ends apart horizontally. Let 
ri (6) and hiO) be two unit vectors attached to the circular strip, where ri points around the circumference, and r2 points along 
the radius. Then a class of low-energy deformations of the paper strip is produced by introducing some bend at every angle 9 
that rotates only around the current direction of n- Denote by u!{9) the rate at which the strip rotates around r2 as a function of 
the angle 9. In a laboratory frame the directions of the unit vectors are determined by 

-^^f2 + uj{e)h (18a) 



h = hy~f2 (18c) 
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FIG. 3; Cut out the circular strip. Rotate the two ends so that the two gray arrows point 180° away from one another, and place them over the 
lower line so that the arrows coincide. The paper will assume a characteristic undulating shape that forms at the edges of leaves or stretched 
plastic, forming a solution close but not identical to that in Fig. ^ with — 2. 



It is easy to check that ri, r2, and r^, remain orthonormal unit vectors under the dynamics described by Eos. |18} This simple 
observation leads to the first important conclusion, which is that the constraint in Eq. and the metric in Eq. |l5|are incompatible 
if the strip undergoes bending alone. To show why, relate coordinate x and angle 9 through 

(19) 



(20) 



/gox = R9 

and note that the location of the center line of the strip is given by 

u{x,0,0) = l{9) = R / d9'h{0') 



However, because of Eq. 18t , one can also write 



l{9)^~Rh{9). 



(21) 



Since r2 is a unit vector, / cannot move more than distance 2R from its starting point. Therefore it is impossible to satisfy the 
constraint in Eq. |l^ when the only bending permitted is around f2- 



VI. SOLUTION WITH CONSTANT BENDING 



Despite the fact that bending about r2 alone cannot solve the physical problem at hand, there is a solution of Eqs. needed 
for further mathematical development. This solution is found when u) is constant. The equations are not obviously solvable even 
in this simple case, since they are nonlinear, but the solution can be expressed in terms of elementary functions. Let 



ri = -Re[— ^ ] (22a) 

r2=A + Re[Be''^/^o], (22b) 
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where A is real but B is complex. With this choice of ri , one immediately satisfies 1 8t . In order for ri to be a unit vector, 

BB = and B-B*=2el. 



Now substitute Eqs. E2 into Eqs. 18c and 18a. They are satisfied if 



B 



Re[iLe'»/»o] ^ ^ e"^/OoBxA-c.c. + 2BxB* +A + Re[Be'''^'^'>] 
9q 46*0 L J 



Matching up coefficients of different powers of exp[i9/0Q] gives 



lUJ - 



B{0,'-1) 



-lUJ - 



290 



-BxA. 



Substituting Eq. 25a into Eq. [25b| and using Eq. Ogives 



(23) 
(24) 

(25a) 
(25b) 



(26) 



Choosing any B that satisfies Eq. ^ one therefore has a solution, where A is determined by Eq. 25a . In particular, when ri (0) — x 
and ^2(0) — y, one has explicitly 



and 



h ^ co^{9/9o)x + 6*0 sin(6'/6'o)y + 9olo 5m{9/9o)z 

h - -0osin(0/0o)i+ (1 - Ol{l - co5{9 / 9Q)))y ~ Lo9l{\ - cos{9/9o))z 

h = -9QUJS.m{9/9Q)x-9luj[\-cos{9/9Q)]y+{9l + [\ - 9^^cos{9 / 9o))z. 

I' ^Rsm{9/9o)9o 

P' ^R{l-cos{9/9o))9l 

l~ ^Roj{l-cos{9/9o))9l 



(27a) 
(27b) 
(27c) 

(28a) 
(28b) 
(28c) 



VII. SOLUTION WITH CONSTANT BENDING AND TORSION 

Since bending around r2 is not a general enough deformation of the strip to satisfy constraint Eq. |T^, one must proceed next to 
consider torsion; that is, bends around r\ . Let describe the rate at which twisting around ^2 occurs, and let ijj describe the rate 
at which twisting around ri occurs. A first case allowing exact solution is when and ijj are constant. Then Eqs. become 

dri 

^cf)rixf2 + h (29a) 

o9 
dr2 

= iPriXr2-ri. (29b) 

o9 



Define 



where 



S2 J LO 4> J \r2 



(30) 



02 + ?/,2 (31) 



Rewriting Eqs. E9^ in terms of these new variables, they become 



dsi ^ ^ ^ 

— =i2+wiiXi2 (32a) 
o9 

= (32b) 

O0 



Eqs. 18 and Eqs. |3^ are identical, except that the latter involve s instead of r. Therefore, choosing a coordinate system where 
x' points along fi (0) and y' points along ^2(0), the i', y' and z' coordinates of si and S2 are given once again by Eqs. ^ 
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VIII. SOLUTION WITH OSCILLATING BENDING AND TORSION 

This solution still does not solve the physical problem at hand, because the strip twists endlessly around the original i axis. In 
order to avoid twisting the strip, the torsion tp must oscillate between positive and negative values. A solution of this type can 
be found by taking 

(j> = ujcosa6 (33a) 
ip = ujsina9 (33b) 

Now the equations describing rotations of the strip unit vectors are nonlinear and have non-co nstan t coefficients. However, 
they can still be solved. Once again, define a new coordinate system with Eq. Substituting Eqs. |33b| into Eqs. |l8] gives now 

—^=S2il+a)+LJSiX S2 (34a) 
^ = -Si{l+a). (34b) 

Defining 

0' = {l+a)0, uj' = Lj/{l+a) and 9'^^ = — (35) 

l+oj' 

one again recovers Eqs. |l8], and its solution Eqs. ^ in terms of primed variables. That is, 

f I = cos{9'/0'o)x + 6*,', sm{9'/0'o)y + O'qJ sin(6l76'^,)f (36a) 

f 2 = -0^ sin ( ^ )i + ( 1 - ( 1 - cos( ^ ) ) )y - ( 1 - cos ( ^ ) )f (36b) 

6^0 y,, 

f3 = -e',J^m(^-)x-e[^uj\\-co^{^-)]y+{e^ + [1 - e[^]co^{^-))z. (36c) 

Finally the solutions can correspond to the shapes at the edge of wrinkled sheets. To see when they are physically acceptable, 
one has to invert all the linear transforms and write the expression for ri explicitly. It is 

f;" = cosa6icos(6'76'^) + 6l^)Sina6'sin(6'76lo) (37a) 

r^' = 6locosa6'sin(6l76'^,) - sina6l(l - 6l(f (1 - cos6l76lo)) (37b) 

=6';;w'cosa6lsin(076'^))+w'e'^'sina6'(l-cos(076'^,)) (37c) 

To satisfy the constraint Eq. |6[ when one integrates f\ to get T, the x component must increase indefinitely, while the y and z 
components must oscillate. The solution acts this way if and only if 

a6 ^ ±e'/e'a (38) 

The two signs in Eq. |8|produce identical solutions, as Eq. |3^ is invariant under O'q — > — 6*0. Adopting the minus sign, one has 

a = -{\+uj^)/2 (39) 
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(40) 



The vectors ri . . . r3 are 



r\ = cos^{a0) - sm^{ae)e'o (41a) 

= [1 - cos(a6l)] sin(Q6')6',f - cos(a6i) sin(a6i)6lo - sin(a6i) (41b) 

[1 - cos(a6l)] sm{ae)9o^ - cos(a6i) sin(a6l)e'^) J (41c) 

r\ = cos(a6') sin(a0)6'o + cos(a6') sin(a6i) (41d) 

r^== {cQ^^{ad)-CQi{aB)\Q'^ -^\x^{ad)B'Q^CQ^(ae) (41e) 

f\ = (^[cos-(a6i)-cos(a6l)]e',f-sin^(a6l)6'^,)cj' (41f) 

r\ = sin(a6l)6'^w' (41g) 

r^ = (cos(a6')-l)6',',V (41h) 

= - cosia0)e'f^ + cos(a6i) + (41i) 
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The wavelength of the pattern is given by the angle through which travels so that the smallest nonzero Fourier component 
of the pattern goes through one period, meaning that 6 travels through In/a. One now finds that the wavelength of the pattern 

A = Z(27r/a) is 



|q;| [l+uj-^y 
Since the strip travels an arclength lirR/a when 9 goes through In/a, one can find L/X= y/go- 



Furthermore, using Eq. nSl one can determine the wavelength A from 



27r 



go 



A 2R{^,-1) 4(V^-1) 



dy 



(42) 



(43) 



(44) 



go and 



Since u; and R, or equivalently A, are the only parameters describing the solution, it has been determined completely by 
the slope of the metric. In particular, the wavelength of the pattern has been determined. 

The functional forms of bending and torsion angles in Eq. 33b were simply pulled out of a hat. One must wonder whether in 
fact the resulting shapes minimize the energy functional in Eq. |10| with which the problem began. This question will be studied 
in later sections, and the answer is "sometimes." That is, there can exist energy functionals and boundary conditions for which 
these solutions are energy minimizers. However, for the particular functional in Eq. |l4| the solutions found in this section are 
excellent approximations but never exact energy minimizers. Furthermore, even for those functionals whose energy they can 
minimize, they do so only for a restricted subset of boundary conditions one can naturally apply to a strip. Returning to the 
paper strip in Fig. ||, the solution described by Eq. ^is only legitimate when the two ends of the strip are placed at a particular 
distance as specified in the caption. However, if one holds the strip in ones hands, it is easy to slide its ends horizontally back and 
forth. The energy- mini mizing shapes resulting from this process certainly exist, but happen not to be in the analytical family 
following from Eq. 33b Finally, one can ask what happens if one takes a strip of length R and allows it to curl up into an 
energy-minimizing shape. Will it ever adopt a periodic structure with the period given exactly by Eq. p4[ ? The answer to this 
question is probably "No," but it has not yet fully been settled 

Some images of solutions from Eq. ^ appear in Fig. The radius R is adjusted in each case so that the wavelength A remains 
constant. Drawing a number of such pictures, one finds that there is an upper limit of « 5.61 for solutions of this type. 
When y/go is larger than this value, the solution collides with itself in the center of the wave and becomes self-intersecting. 




FIG. 4: Height of edge /' (0) versus length of edge / '^(S) for solution given in Eq. ^ and various values of 



go- 



IX. IMPLICATIONS FOR ENERGY FUNCTIONAL 

Assuming that displacements are of the form given in Eq. |3^, it is possible to return to the functional Eq. [l^ and decide which 
terms dominate the energy of a stretched sheet. From the customary theory of thin sheets one expects the energy to come from 
the sum of two terms. The first term is the stretching energy of the sheet, and is of order f , where t is the thickness of the sheet; in 
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Eq. |l^ the stretching energy would arise from the terms on the top of the right hand side that involve no derivatives with respect 
to z- The second term is the bending energy of the sheet, and it is proportional to P. 

This expectation is correct. Bending and stretching energies dominate the energy functional. Which of them dominates is a 
subtle matter The stretching energy is proportional to the thickness of the sheet f while the bending energy is proportional to 
the thickness of the sheet cubed, t^. In the limit of thin sheets, one would think that the stretching energy would be dominant. 
However, matters are not that simple, because the energy is proportional to the width of the strip w to the fifth power, . The 
limit 1 believe is most interesting in practice is 



t W 



(45) 



In this limit the bending energy proportional to t^/R^ dominates the energy, while the stretching energy, proportional to w*/R'^, 
can be neglected. 

To obtain the energy formally, take u in the form 



0)^2 - 






h/R 


(2)^2 _ 






h/R 






^q^f] 


h/R 



(46) 
(47) 
(48) 
(49) 



with related to x by Eq. |l9|. 

If all the coefficients q were zero, one would have 



du 
dx 



go d 



R 86 



l{0) + hy + hz 

y dh z dh 
R 00 R 86^ 



du 
dy 
du 

Tz^''- 



(50a) 

(50b) 
(50c) 

(50d) 
(50e) 



However, the functions q are needed to allow small contractions of the membrane so as to minimize its energy. The minimiza- 
tions are performed in the following way; first find the term in Eq. [l^ proportional to y^ /R^- This term can be made to vanish 
through a unique choice of q^zy , q^zy , qfy , qyy ^and qfy , where q^^y — qzf — ~'^qyy , and the others vanish. The next term to 
consider is the term proportional to z^ /R^, as this is larger than the term proportional to y'^ /R^- Minimizing the term proportional 
to determines the following values for the variables q: 



q^ 




q^ 






q^ 


= V2(/)/6, 


qS 


= 0/3, 


?g' = o 


qS 


= 0/12, 


q^ 


= 0, 


?(:;' = -0/6 



(51a) 
(51b) 
(51c) 

Inserting the expressions for into Eq. |o[ one has the following results. Since there are are only two layers in the z 
direction, carry out that sum explicitly, but convert sums in the x and y direction to integrals, using the fact that the area per 
particle is ^J2>/Aa^. The leading term proportional to the sheet thickness is of the form 

K — T / dx [many terms!] (52) 

R L Jo 

but the term that dominates the energy when inequaUties Eq. ^ hold is 

at^Lwl f\ (W^ + 2gg0^) 
/C— 5— - / dx- ;= -. (53) 



R^ LJo V3 
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When Eq. ^ dominates the energy, the strip is wide enough so that t, but not so wide that it is favorable to begin buckling 
or forming fine structure in the y direction. 

The various powers of go appearing in Eq. ^ arise because the microscopic lattice underlying the calculation has been 
stretched, and its elastic properties made anisotropic by the new metric. The calculations are producing a specific description for 
how stretching material in the x direction by a factor of ^/go changes its elastic properties. This description is unlikely to apply 
to the polymeric materials in which most experiments are conducted. For this reason, direct numerical studies of the bending of 



strips have been carried out in such a way that go in Eq. 53 should be set to 1. Rather than starting with a strip of length L—\, 
stretching it to length X^/go, and then introducing a gradient in the metric along y, the numerical studies start with a strip of 
length L, introduce a linear gradient of slope l/Rin the metric along y, and then constrain the ends of the strip to sit at distance 
A = L/^/go- That is, the constraint at the ends of the strip instead of Eq. 16 is 



The energy of a strip prepared in this way is 



ii{x + L,y,z) —xX + u{x,y,z)- (54) 



LJo V3 

X. EULER-LAGRANGE EQUATIONS FOR MINIMUM OF ELASTIC ENERGY 



The solution of Eq. |3^ was obtained by assuming that <j) and ijj have the forms given in Eq. 33b These choices for the torsion 
and bending of the strip did make it possible to find a shape for the strip that was consistent with the constraint in Eq. [l^, and 
with a few other physical considerations. 

Varying (j) and ip away from the arbitrary forms obtained by guessing leads to additional solutions, generally of lower energy. 
Depending upon the precise form of the continuum energy functional, and the boundary conditions to which the problem is 
subjected, the solutions in Eq. ^ can turn out to be exact, and in many other cases are excellent approximations to exact 
solutions. 

To look for energy-minimizing shapes, 1 will use Eq. |5]to express the energy of a strip, and will rewrite it as 

£ = Jde^^^ + ^i;\ (56) 

If Ci is equal to C2 then the analytic solution of Eq. ^can minimize the energy for certain boundary conditions. According to 
Eq. |5[ in fact the constants are in the ratio of 2 to 3, and since they are not equal, the solution is approximate, although for some 
boundary conditions the approximation is excellent. 

In explaining the calculations leading to these conclusions, it is useful to derive the equations of motion Eq. |l8| in a different 
way. Let ri (0) ... h [6) be orthonormal vectors that evolve as a function of 6. Letting primes denote derivatives with respect to 
0, one has in general that 

r' = Yyrrj)n (57) 

./ 

This equation of motion is subject to a number of constraints. First, each vector must retain unit magnitude. Therefore 

|r,|2 = 0^r,'-r,- = 0. (58a) 

Second, since the vectors are orthogonal, for / ^ j 

r, • fj = 0. (58b) 

Third, the condition that ri and r2 curl around one another because of the gradient of metric, but otherwise be rigid against 
twisting around means that 

ri-r^ = -l or r[-r2 = \. (58c) 

There are no constraints upon h ■ r[ or ^3 -^2; define these functions to be 

^ = h-f[ (59a) 
ijj = h-r'2. (59b) 
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Then Eqs. |l^ follow from Eq. ^ just by employing Eqs. p8|and |59| . 

Once vectors rj . . . ^3 are subject to the constraints in Eqs. |58|, the energy Eq. 
constant, as 



Jde^\r[\' + ^\r',\' 



)qcan be rewritten, apart from an overall additive 



(60) 



Imposing the final constraint 



one needs to extremize the functional 



/' 

Jo 



dB^\r[\^ + ^\r',\^ + p-h(e)-\umr]\- «2(^?)n 



u,{e)\rl\-u,{e)r, 



(61) 



(62) 



The constant vector p and the functions ui{9) . . . u^{9) are Lagrange multipliers enforcing the constraints. 

There may be some concern over the legitimacy of using equations of constraint in order to simplify the form of the energy 
functional before carrying out the variation. The results will not be reported here in detail, but the entire variational procedure 
has been repeated using ?3 • r[ and r-^ ■ r'2 as the definitions of (p and ip, and enforcing directly the constraints r\ = (pr^ xr2 + ^2 
and rj = x ''2 ^ n ■ The algebra is more involved than what is reported below, but the results are identical. 

Taking the variation of the functional Eq. ^ with respect to ri and r2 gives 

— C\r\' + p — uiTi — U2r2 — U4r2 — 
-Car" - "2?! - U3r2 + U4r[ + u\ri = 0. 

Now employing Eq. |l8[ one can rewrite Eqs. |63b| as 

Ci [-(0^ + l)ri - (pxpfi + [4)' + V']] + u\r\ + U2r2 + M4[V''^3 - n] = /? 

C2[-(t)'^f\ - [il? + 1)^2 + W - 0)^3] + M2n + 1(3^2 - M4[0r3 + ^2] - u'^Tx = 



(63a) 
(63b) 



(64a) 
(64b) 



One obtains six scalar equations by taking the dot product of Eqs. 64 with f\, r2 and ^3 in turn, using the fact that they are 



orthonormal once all the Lagrange multipliers are chosen properly. The only time that ui appears is when one takes the dot 
product of Eq. 54a with ri, and similarly the only time that M3 appears is when one takes the dot product of Eq. |64b| with ^2- 
Therefore, ui and M3 can be taken to be whatever is needed to satisfy these two equations, which need not even be written down. 
The remaining four equations are 



U2^C2(j>1p + U^ 
U4(f) — C2[tp' — (p] 

U2^p-r2 + Ci4'ip 

U41I} ~ p ■ f3~C][ij)' + 1p] 



(65a) 
(65b) 
(65c) 
(65d) 



Eliminating U2 allows one finally to write a complete set of Euler-Lagrange equations determinining minima of the energy: 

u\=p-r2 + [Ci- C2] (pip (66a) 

ip' = U4ip/C2 + 'p (66b) 

^' =[p- {ri 'xh)- U41P] /Ci - ^ (66c) 

r'l = (/)(ri X ^2) + h (66d) 

r2 = -ijj{ri X ^2) - h . (66e) 

Notice that all the terms in Eq. ^have definite parity under 0^—0. Assuming that all energy minima have the same symmetry, 
p must point along x, (p> and M4 must be even, and t/j must be odd. The solutions of Eq. ^ are indexed by four constants, the 
initial values of (f) and M4, and by p ^, and finally the angle OfXo which the solutions are to be integrated. 

To thi s point, the manipulations make use o f Eq . but do not use any features of the particular forms of (p ™d i/) chosen in 
Eq. 33b. Do (p> and t/j chosen according to Eq. [33b| solve Eqs. ^? Inserting these forms of and ^, one quickly finds that 

(67) 



M4 = C2(a— 1), and u^ = 

(C2 — C\)uj^ sin a9 cos aO — p-r2 
(Ci — C2)(l — Q;)wsin q;6' = p-h- 



(68) 
(69) 
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(a) 



Distance x 




FIG. 5: (a) Comparison of three separate methods to find shape of strip. The shape indicated as the exact solution results from solving Eq. |66| 
with M4(0) = -9.26053, (/>(0) = 1.8129, and px = -4.4599, leading to a solution with ^/g^ = 1.333 and Of = 3.1416. The shape indicated 
as an approximate solution comes from evaluating Eq. ^ for ^/gQ = 4/3, 9f = tt, and R = 63.7 = 200/7r. Finally, the shape resulting from 
direct simulation is obtained as a result of a direct numerical minimization of Eq. MOl for a strip 201 atoms long in the x direction, and 12 layers 
wide in the y direction, constrained to have horizontal period of 150 so that y/go = 4/3, and R = 63.7. The energy predicted by the function 
Eq. ^ with L — 200 and w = II ^3/4, and the approximate forms for (j> and ip from Eq. 33b is 1 .357 1 A^^a"' . The corresponding calculation 
using exact results from Euler-Lagrange equations in Eq. pq is 1.3557K;a^ The energy found from direct numerical minimization of Eq. [KJ is 
1 .377 ICa^ . (b) Three-dimensional visualization of the final shape resulting from the numerical minimization. 



Turning to Eq. ^l| one finds that p must point only in the i direction and 

{C2-Ci)u;^=Px{l+e'o) 
(C,-C2)(l-a)c^ = M)^'- 



Using Eq. B3 to relate 6^ and a, and Eq. B5l one finds 



Ci =C2 



(70a) 
(70b) 



(71) 



a = 1 + 



(72) 



Unfortunately, Eq. ^ contradicts Eq. Therefore, the functional forms in Eq. 33b minimize the energy only if Ci = C2, and 
otherwise do not. 



XI. NUMERICAL SOLUTIONS 



The Euler-Lagrange equations in Eq. 66 constitute a closed set of equations that determine all energy minima. It is surprisingly 
tricky to employ them for this purpose, and 1 have only partial results to report. Searching for periodic solutions of Eqs. ^ 
proceeds in the following way: choose a value of p^, and the initial values of U4{9) and 0(6'). Integrate forward in 6 until 
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Distance x 



FIG. 6: Comparison of the lowest-energy shape from Fig. |g with the lowest-energy shape found numerically by varying the initial conditions 
in Eqs. ^ while requiring L/ A = y/gQ = 4/3 and scaling the solutions by i? = 200/7r. The longer-wavelength solution has energy 1 .2075/Ca^ . 
Lower-energy solutions probably exist, but my numerical routines have been unable to track beyond this point. 



reaching an angle 0f at which ip{9f) vanishes. It will not in general be true that the solution at this point obeys Eq. |l6[ or that 
ri and r2 point along the i and y axes as they must. However, by adjusting p^^ alone, all of these boundary conditions can be 
achieved at once. One still has the two initial values </)(0) and M4(0) to vary. These allow one to vary 9/ and A independently. 

I have used a shooting algorithm to find solutions along these lines. The process works, but is very delicate. Unless all 
constants are chosen very near to their final values, the solutions of the equations fly off into completely unphysical regions. 
The analytical solutions of Section VIII are indispensable for the purpose of locating good starting points for shooting solutions. 
Once a few solutions have been located, additional solutions can be obtained by using extrapolation to predict new values for 
initial conditions and searching in a narrow range of parameter values around the extrapolated predictions. 

When C\ — C2 in Eg. |56|, the solution in Eq. ^ does indeed minimize the energy when the wavelength of the pattern is 
constrained to obey Eq. 44]. However, even when Ci = C2, the lowest-energy solution for a given is not given by demanding 
that the wavelength correspond to this analytically tractable value. For a given ^/go = L/\, energies per length around 25% 
lower are obtained by allowing A and L to increase; the energy decreases monotonically as they become larger. The shooting 
procedure eventually becomes unable to track additional solutions, and it is not yet possible to say whether there are any energy 
minima to be found at particular wavelengths. 

When Ci =2 and C2 = 3 in Eq. as it should for an isotropic two-dimensional material in accord with Eq. the solution 
in Eq. ^ provides a good enough approximation to the true solution that the numerical routines are able to converge by using 
its properties as a starting point. The comments of the paragraph above otherwise do not need to be changed. As shown in Fig. 
|[ when boundary conditions corresponding to a solution obeying Eq. are imposed upon the Euler-Lagrange equations, their 
solution is almost indistinguishable from Eqs. However for the same value of y/go, the longer the value of the wavelength 
A, the lower the energy. There is no indication of an energy minimum at a periodic solution. Fig. ^compares the lowest-energy 
solution I have found for = 4/3 with the solution displayed previously in Fig. ||. 

An alternate way to check whether periodic solutions are favorable is to carry out direct numerical minimization in large 
systems. Returning to the numerical setting of Fig. ^ I triple the length of the system, fix the boundaries at distance 3A, leave 
all else as before, and minimize the energy Eq. [T^. The strip develops a composite structure in which the long wavelength is 
3 times the short wavelength, and with an energy that is somewhat less than 3 times the energy of the structure 1/3 as long, as 
shown in Fig. ^ 

Thus it is possible that the lowest-energy state for the simple metric in Eq. |l5| involves a cascade of oscillations on many 
scales, but the matter is not yet settled. Many other problems remain to be addressed, including the application of the ideas 
obtained here to the more complex metrics illustrated in Fig. 
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Distance x 

FIG. 7: Numerical minimization of the same system described in Fig. ^ but with A and L three times as long. The total energy is 

which is less than the value 3 x X.'ilK.cr' one would obtain simply by pasting together three copies of the numerical solution in Fig. ||. This 

computation provides an additional hint that periodic solutions do not minimize the energy. 



Acknowledgments 

Eran Sharon discussed these problems with me for weeks, and spent months carrying out experiments, before I finally began to 
see there was something interesting to be done. Lorenzo Sadun helped us see what could be gained from differential geometry. 
Thanks to Ralf Stephan for pointing out Ref. [|]]. Financial support from the National Science Foundation (DMR-9877044 
DMR-0101030) is gratefully acknowledged. 



1. E. Sharon, B. Roman, M. Harder, G.-S. Shin, and H. L. Swinney, Nature 419, 579 (2002). 

2. M. Byrne, M. Timmermans, C. Kidner, and R. Martinssen, Current Opinion in Plant Biology 4(1), 38 (2001). 

3. R Green, C. Steele, and S. C. Rennich, Annals of Botany 77, 515 (1996). 

4. R B. Green, American Journal of Botany 86, 1059 (1999). 

5. S. Nechaev and R.Voituriez, Journal of Physics A 34, 11069 (2001). 

6. D. Henderson and D. Taimina, Mathematical Intelligencer 23, 17 (2002). 

7. A. V. Pogorelov, Differential Geometry (P Noordhoff N. V., Groningen, 1956). 

8. L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Pergamon Press, Oxford, 1986), 3rd ed. 

9. E. H. Mansfield, The Bending and Stretching of Plates (Pergamon, New York, 1964). 

10. H. S. Seung and D. R. Nelson, Physical Review A 38, 1005 (1988). 

11. A. E. Lobkovsky and T. A. Witten, Physical Review E 55, 1577 (1997). 



